\documentclass[10pt]{article}
\usepackage[hyperref,UTF8,scheme=plain]{ctex}
\usepackage[left=5.00cm, right=3.00cm, top=5.00cm, bottom=5.00cm]{geometry}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{float} % 更好地控制图片位置
\usepackage{listings}
\usepackage{xcolor}
\usepackage{framed} % 可以用 \begin{shaded}，即背景色块
\definecolor{shadecolor}{rgb}{0.9,0.9,0.9}

\hypersetup{
    colorlinks = true,
    linkcolor= blue
}

\lstset{
	basicstyle=\footnotesize,
	numbers=left,
	breaklines,%自动换行
	frame=shadowbox,%加边框
	backgroundcolor = \color{red!30!green!30!blue!30},
	%背景色：蓝灰色
	language=C++, % 语言为 c++
	commentstyle=\color{red!80!green!80!blue!80},% 注释用蓝灰色字体标出
	keywordstyle=\color{blue!70} % 高亮部分用蓝色字体
}

\begin{document}

{\Huge \bfseries 2019年秋《计算物理》期末试题答案}\\

\vspace{1cm}

\section{单项选择题，每题3分，共10题，总分30分。}

\kaishu

ADDDD BDDDD

\section{解答题，每题10分，共7题，总分70分}
2.1 请写一个gnuplot脚本，制作绳子上驻波的动图：
\begin{equation}
u(x,t) = \cos(2x) \cos(0.1*t),
\end{equation}
其中$u(x,t)$表示$x$处质点在时间$t$的位移，使用国际单位制。要求：
\begin{itemize}
\item[a.]动图纵向坐标范围是$[-1,1]$。
\item[b.]横向坐标范围是$[0,10]$。
\item[c.]时间范围是$0\sim 80$秒。
\end{itemize}
答案：
\begin{lstlisting}
set term gif animate
set output "2.1.gif"
set xrange [0:10]
set yrange [-1:1]
do for[t=0:80]{
        plot cos(2*x)*cos(0.1*t) w l lw 3
}
set output
\end{lstlisting}
{ \color{red} 评分标准：成功生成动图(7分)，横坐标、纵坐标范围(2分)，时间范围(1分)。}
\vspace{0.3cm}

\noindent
2.2 请编写1维Ising模型代码，用周期性边界条件、100个原子模拟1维晶格。其哈密顿量为
\begin{equation}
H = - J \sum_{\{i,j\}} s^{(i)}_z s^{(j)}_z,
\end{equation}
其中$J$表征相邻原子相互作用，可以设为能量单位，即设为$1$，$\{i,j\}$表示相邻的两个原子。
请用Markov链+Metropolis模拟，计算平均每原子自旋第3分量$\langle s_z \rangle$与$kT/J$的关系，作图并给出结论：1维Ising模型存在相变吗？如果有，相变点为$kT/J = $？

答案：不存在相变。代码如下：
\begin{lstlisting}
/*
 * 1-dimensional Ising model, with Markov chain in Monte Carlo method
 *      H = - J \sum_{ij} sz_i sz_j - muB \sum_i sz_i,
 * where J is the interaction strength between two atoms, negative/positive for ferromagnetic/anti-ferromagnetic material, k is the Boltzman constant, T is the temperature, we set J as unit.
 */

#include<iostream>
#include<stdlib.h>
#include<time.h>
#include<fstream>
#include<cmath>
using namespace std;

void Ising1d(int length, int L, double kT, double muB, double & sz_ave, double & e_ave, double & e2_ave){

        long i,j,k, Sz=0, Sz_neighbour, count=0;
        double J=1, E=0, dE, sz_sum=0, e_sum=0, e2_sum=0;
        int *spin=new int [L];//spin of each atom, which can only be 1 or -1
        bool flag_flip;

        srand((unsigned)time(0));//initialize rand() with a random seed, according to the real time when this code runs.
        for(i=0;i<L;i++){
                spin[i] = 1 - 2 * (rand() %2);
                Sz += spin[i]; // initial total spin
        }
        // calculate initial energy
        for(i=0;i<L;i++){
                if(i==0) Sz_neighbour = spin[1] + spin[L-1];
                else if(i==L-1) Sz_neighbour = spin[L-2] + spin[0];
                else Sz_neighbour = spin[i-1] + spin[i+1];
                E += - J * spin[i] * Sz_neighbour - 2 * spin[i] * muB;
        }
        E /= 2;

        for(j=0;j<length;j++){

                for(i=0;i<L;i++){

                        Sz_neighbour = spin[ i==0 ? L-1 : i-1 ] + spin[ (i == L-1)? 0 : i+1 ];

                        dE= 2 * J * Sz_neighbour * spin[i] + 2 * muB * spin[i];

                        flag_flip=false;
                        if(dE<0){
                               flag_flip=true;
                        }
                        else{
                                if( (double)rand()/RAND_MAX < exp(-dE/kT) ){// rand()/RAND_MAX uniformly distributes in [0,1], so it has possibility exp(-dE/kT) to be less than exp(-dE/kT)
                                        flag_flip=true;
                                }
                        }

                        if(flag_flip){
                                E += dE;
                                spin[i] *= -1;//flip the spin
                                Sz += 2*spin[i];// total spin and total energy changes
                        }
                }
                // after one round, collects data

                if(j>0.2*length){// We assume the system is in equilibrium after 0.2*length

                        /*
                        for(i=0;i<L;i++){
                                if(spin[i]==1)cout<<"+";
                                else cout<<"-";
                        }
                        cout<<"\tSz="<<Sz<<"\tE="<<E;
                        cout<<endl;
                        */

                        sz_sum += (double)Sz/(L);
                        e_sum += E/(L);
                        e2_sum += E*E/(L);

                        count ++;
                }
        }

        sz_ave = sz_sum/count;
        e_ave = e_sum/count;
        e2_ave = e2_sum/count;

        delete [] spin;
}

int main(){

        int length = 1E4, L=100;
        double sz_ave, sz_ave_sample, e_ave, e_ave_sample, e2_ave, muB=0;

        ofstream fp("Ising1d.txt");
        for(double kT=0.1; kT<=1; kT+=0.02){
                sz_ave = 0; e_ave = 0;
                for(int i=0;i<10;i++){
                        Ising1d(length, L, kT, muB, sz_ave_sample, e_ave_sample, e2_ave);
                        sz_ave += sz_ave_sample;
                        e_ave += e_ave_sample;
                }
                sz_ave /= 10;
                e_ave /= 10;
                fp<< kT <<"\t"<< fabs(sz_ave) << "\t" << e_ave <<endl;
        }
        fp.close();

        return 0;
}
\end{lstlisting}
\begin{figure}
\centering
\includegraphics[width=0.8\linewidth]{../src/final/2}
\caption{一维Ising模型，平均磁化率与$kT/J$的关系。}
\label{fig:2}
\end{figure}
{\color{red} 评分标准：写出完整的1维Ising模型代码(5分)，正确处理马尔科夫链($dE \rightarrow P(X -> X')$)(3分)，用图片展示计算结果(2分)。}
\vspace{0.3cm}

\noindent
2.3 已知两个物理量$X, Y$之间的关系为$Y = k X + b$，现有十组实验数据
$(X,Y) = (1,3.2), (2,5.1), (3,6.3), (4, 7.01), (5, 7.05), (6, 8.5), (7, 9.2), (8, 10.8), (9, 12.1), (10, 13)$，
请用最小二乘法拟合出$k, b$的值。

答案：运行如下代码，得到 $ k = 1.01648,	 b = 2.63533$，相应地，$RMSD=0.425373$。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

#include"library.h"

double g0(int n, double *x){
        return 1;
}

double g1(int n, double *x){
        double z = 0;
        for(int i=0;i<n;i++) z += x[i];
        return z;
}

double g2(int n, double *x){
        double z = 0;
        for(int i=0;i<n;i++) z += x[i] * x[i];
        return z;
}

double g3(int n, double *x){
        double z = 0;
        for(int i=0;i<n;i++) z += pow(x[i],3);
        return z;
}

int main(){

        int i,j,n=1, m=2,N=10;
        double ( *g[3] )(int, double * )={g0,g1,g2};

        double **x_data = new double * [N];
        for(i=0;i<N;i++) x_data[i] = new double [n];
        double *y_data = new double [N];

        double DATA[10][2] = {{1,3.2},{2,5.1},{3,6.3},{4,7.01},{5,7.05},
        {6,8.5},{7,9.2},{8,10.8},{9,12.1},{10,13}};

        for(i=0;i<N;i++){
                x_data[i][0]=DATA[i][0];
                y_data[i]=DATA[i][1];
        }

        double *a = new double [m];
        cout<<"RMSD="<<LeastSquare(N, x_data, y_data, n, m, a, g, 0.01 )<<endl;
        cout<<"a: "; for(i=0;i<m;i++)cout<<a[i]<<","; cout<<endl;

        delete [] a, y_data;
        for(i=0;i<N;i++)delete [] x_data[i]; delete [] x_data;

        return 0;
}
\end{lstlisting}
{\color{red} 评分标准：写出完整的最小二乘法代码(5分)，正确选取理论公式(3分)，输出正确计算结果(2分)。 }
\vspace{0.3cm}

\noindent
2.4 请任选一种数值方法，求解以下势井的基态能量，
\begin{equation}
V(x) = \left\{
\begin{aligned}
& x^2, ~~ 0<x<1, \\
& \infty, ~~ else.
\end{aligned}
\right.
\end{equation}
我们设定$\hbar = m = 1$，体系定态薛定谔方程为
\begin{equation}
- \frac{d^2}{dx^2} \psi(x) + V(x) \psi(x) = E \psi(x).
\end{equation}

答案：使用龙哥库塔方法+打靶法，运行下面的代码，得到基态能量为$E_{g.s.}=5.21525$。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>
#include"library.h"

double E;

void derivative(int n, double x, double *y, double *dydx){
        dydx[0] = y[1];
        dydx[1] = 2*(x*x - E)* y[0];
}

double psib(double x){

        E=x;

        int N=10000,n=2;
        double a=0, ya=0, deriv_a=1, b=1;
        int i,j;
        double h = (b-a)/N;
        double *y = new double [n];
        double *ynext = new double [n];
        y[0]=ya; y[1] = deriv_a;
        for(i=0;i<N;i++){
                rk4(a+i*h, h, n, y, ynext, derivative);
                for(j=0;j<n;j++) y[j] = ynext[j];
        }
        return y[0];
}

int main(){

        double eig = bisection(psib, 0, 20, 1E-9);
        cout<<"E="<<eig<<endl;
        // g.s. E = 5.21525

        return 0;
}
\end{lstlisting}
{\color{red} 评分标准：写出完整的最小二乘法代码(5分)，正确选取理论公式(3分)，输出正确计算结果(2分)。 }

\noindent
2.5 请产生服从如下分布的离散随机数
\begin{equation}
P(x) = \left\{
\begin{aligned}
&0.5, ~~ x = 1,\\
&0.5, ~~ x = -1.
\end{aligned}
\right.
\end{equation}
进行10000次抽样，得到它们的平均值
\begin{equation}
\overline{x} = \frac{1}{10000} \sum^{10000}_{i=1} x_i,
\end{equation}
试验证：$\overline{x}$值大致服从正态分布$N(\mu = 0, \sigma^2 = 10^{-4})$。\\

答案：运行如下代码，作图验证中心极限定理
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cstdlib>
#include<cmath>
#include<fstream>

#include"library.h"

double sample(int n, double a, double b, double (*func)(double x), double fmax){

        double y=0;
        for(int i=0;i<n;i++){
                y += (rand()%2) * 2 -1 ;
        }
        return y/n;
}

/*
 * count_random(...) 统计随机数rand出现在 [a,a+h), [a+h, a+2h), ... 中的频率，记录在count中
 */
void count_random(double a, double h, int * count, double rand){

        count[ (int)((rand-a)/h) ] ++;
}

double f(double x){//[-3:3]之间均匀分布
        if(x >= -3 && x <= 3) return 1.0/6;
        else return 0;
}

double NormDistr(double mu, double sigma, double x){
        return 1 / sqrt(2*M_PI) / sigma * exp( - (x-mu)*(x-mu)/2/sigma/sigma );
}

int main(){

        int i,N=10000;
        double a=-0.1, b=0.1, fmax = 1.0/6;
        double h=0.001, x;
        int n = (b-a)/h;
        int *count = new int [n];
        for(i=0;i<n;i++) count[i] = 0;
        for(i=0;i<N;i++){
                x = sample(10000, a, b, f, fmax);
                count_random(a, h, count, x);
        }
        ofstream fp("2.5.txt");
        for(i=0;i<n;i++){
                fp<< a+ (i+0.5)*h <<"\t"<<(double)count[i]/N/h<<"\t"<< NormDistr(0, 0.01, a+ (i+0.5)*h )<<endl;
        }
        fp.close();
        delete [] count;

        return 0;
}
\end{lstlisting}
\begin{figure}
\centering
\includegraphics[width=0.7\linewidth]{../src/final/5}
\caption{中心极限定理的验证：如2.5题设中的分布取随机数，取10000次抽样，抽样平均服从正态分布。}
\label{fig:2.5}
\end{figure}

\noindent
2.6 请任选一种方法，完成如下数值积分，并与解析解比较
\begin{equation}
I_f = \int^1_0 e^{2x} dx.
\end{equation}

答案：使用梯形法积分，运行如下代码，得到：$I_f \approx 3.19453$。
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

#include"library.h"

double f(double x){
        return exp(2*x);
}

int main(){

        double a=0,b=1;

        cout<<"integ = "<<integ_trapzoid(a,b,1000,f)<<endl;
        cout<<"analytic: If = "<< 0.5*(M_E*M_E-1) <<endl;
        return 0;
}
\end{lstlisting}
{\color{red} 评分标准：声明使用的积分方法(2分)，写出正确代码(5分)，得到正确结果(2分)，与精确值比较(1分)。 }

\noindent
2.7 尝试用VMC方法，用试探波函数$\psi(x) = x e^{-\alpha x^2}$，求解1维势井基态：
\begin{equation}
V(x) = \left\{
\begin{aligned}
&x^2/2, ~~~ 0<x<\infty\\
&\infty, ~~ x<0
\end{aligned}
\right.
\end{equation}
薛定谔方程变为
\begin{equation}
(-　\frac{1}{2} \frac{d^2}{dx^2} + V(x) ) \psi = E \psi.
\end{equation}
理论解为
\begin{eqnarray}
E_{g.s.} &=& 1.5, \\
\psi_{g.s.} &\propto& x e^{-x^2/2}.
\end{eqnarray}

答案：运行如下代码，得到$E^{vmc}_{g.s.} = 1.5$。
{\color{red} 评分标准：写出正确的局域能量公式(2分)，写出完整正确的变分蒙特卡洛方法(5分)，得到正确结果(2分)，文本/代码具有良好可读性(1分)。 }
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<fstream>
#include<cmath>

#include"library.h"

double psi(double alpha, double x){
        return x*exp(-alpha*x*x);
}

/*
 * EL(alpha, x) = (- 1/2 *d^2/dx^2 + x^2/2) psi / psi = 3 alpha + (1/2-2*alpha*alpha)*x*x
 */
double EL(double alpha, double x){
        return 3*alpha + (0.5-2*alpha*alpha)*x*x;
}

int main(){

        ofstream fp("2.7.txt");
        for(double alpha=0.4; alpha <=0.6; alpha += 0.01){
                fp << alpha << "\t" << vmc1d_aveEL(alpha, 100, psi, EL, 0, 5, 0.2, 10000) <<endl;
        }
        fp.close();
        return 0;
}
\end{lstlisting}
\begin{figure}
\centering
\includegraphics[width=0.7\linewidth]{../src/final/7}
\caption{2.7 One-dimensional variational monte carlo method: find the ground state energy}
\label{fig:7}
\end{figure}

\end{document}